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Accurate and reproducible measurements of current- 
voltage characteristics have been reported for single 
molecules i 2 - The experiments have revealed a wealth of 
new interesting physical phenomena in molecules such 
as Coulomb blockade, Kondo effects, negative differen- 
tial resistance, vibronic effects and local heating, as well 
as switching and hysteresis^ Much of the theoretical 
and computational studies have been based on Keldysh 
nonequilibrium Greens functions,— ifiiSui but the under- 
standing of fundamental mechanisms of quantum trans- 
port in nanoscale molecular systems also requires the de- 
velopment of new theoretical approaches to nonequilib- 
rium interacting many-particle quantum systems. In this 
paper, we develop the transport theory, which is based 
on some more recent advances in quantum field theory of 
nonequilibrium states. 

To carry out our program we will use the formalism 
of thermo field dynamics (TFD)^ The TFD is one of 
the approaches to quantum field theory at finite temper- 
ature. It is an operator based formalism which is compli- 
mentary to Matsubara and nonequilibrium Green's func- 
tions. TFD has been successfully applied to nuclear 
physics ) 11 ' 12 ^ 3 particle physics^ and problems of sta- 
tistical mechanics, 15 although its applications to chem- 
ical physics have been very limited. As we discuss in 
this paper, TFD can also be naturally used to describe 
nonequilibrium steady states of many-particle quantum 
systems, in general, and quantum transport of electrons 
through the interacting region, in particular. 

To begin, we first summarize relevant for us aspects 
of TFD. A nonequilibrium quantum system can be de- 
scribed in terms of density matrix p(t) that is obtained 
as a solution of Liouville equation (h = 1): 



ip{t) = [H,p{t)] 



(1) 



We can introduce the square root of the density matrix 
(p(t) — fit)^ (t)), which satisfies the same time evolu- 
tion equation as the density matrix: 



n(t) = [H,j(t)] 



(2) 



The average value of an operator at particular moment t 
can be computed as 



The TFD defines a time-dependent wavefunction \ip(t)) 
in such a way that the matrix element over this wave- 
funcion (if>(t)\A\ip(t)} agrees with nonequilibrium aver- 
age ([3]). We assume here for simplicity that the eigenval- 
ues of the system Hamiltonian H are discrete: H\n) — 
-^n|^); {m\ri) = S nm . TFD introduces a fictitious sys- 
tem (denoted by tilde) which is identical copy of the 
original system: H\n) — E n \h), (rh\h) — 5 nm . Then 
we can write wavefunction for nonequilibrium state as 
|"0(*)) = Y,nrnlrnn{t)\n) ® |m) , where 7„ m = (n|7|m) are 
matrix elements of the square root of the density matrix 
computed in the basis of eigenvectors of the Hamiltonian 
of the system. The average of operators defined in our 
physical space A = A® I over this wavefunction gives 
exactly nonequilibrium average ©. The Liouville equa- 
tion ([2]) is equivalent to the time-dependent Schrodinger 
equation 



im))=n\m), 



(4) 



where H = H — H is the so-called thermal Hamiltonian 
of the system^ TFD should not be confused with sim- 
ple notational reformulation of the Liouville equation. 
There are profound physical differences between the two 
approaches, which will be discussed in the end of the next 
paragraph. 

Let us consider the quantum system which is brought 
to nonequilibrium by interaction with one or several 
macroscopic thermal baths or particle reservoirs. The to- 
tal Hamiltonian H = Hs+Hb+Hsb includes the system 
Hamiltonian Hs , the Hamiltonian of the baths Hb , and 
the system-bath interaction Hs b and it is assumed to be 
time-independent. The corresponding thermal Hamilto- 



nian is 



H = H - H = H S +H B + H sb . 



(5) 



(A{t))=Tr(p(t)A)=Tr(^{t)A 1 {t)) 



(3) 



Here we denote all thermal baths collectively by B and 
all operators which are defined in the bath Hilbert space 
will be also denoted by B, for example, Ab- The equa- 
tion (|3]) should be solved subject to the special boundary 
conditions for operators acting in Hilbert subspace of the 
baths. Namely, at t — > — oo, when the coupling between 
baths and the system Hsb = 0, the baths were prepared 
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in the canonical ensembles at temperatures T B : 

um (m\Mm) = ^-Y, e ~ Eb/TB( t>\ AB \i>)i (6) 

where Z B is the partition function of bath B and \b) are 
eigenstates with eigenergies E\, of the Hamiltonian of the 
bath B. Let us consider the particular case, which will be 
very important for our development of transport theory, 
namely, we consider the baths containing non-interacting 
fermions Hb = J^B^beB £ b a l a b- Here a],(ab) creates 
(destroys) an electron in the state b of the bath B. As was 
demonstrated by Umezawa^, in this case \ip(t — > -co)) 
becomes a vacuum for operators /3b, /3b, which are related 
to the original bath particle creation and annihilation op- 
erators via thermal canonical Bogoliubov transformation: 

ft = vTfn^ab - ) :beB (7) 

Pb = v 1 - n b a b + ^/nt, aj J 

Here rib = (l + eS Eb ~ fJ/B ^ TB )~ 1 is the Fermi-Dirac occupa- 
tion number for level b in the bath and the canonical Bo- 
goliubov transformation should be performed separately 
for all baths B. Now we are ready to discuss the impor- 
tant advantages of the TFD approach over the Liouville 
equation. TFD operates in the space of the distribu- 
tion amplitudes (Hilbert space), whereas the Liouville 
operator is defined in the space of distribution functions. 
As a consequence, unitary representations of states and 
the definition of quasiparticles as unitary transformation 
of the original particles are not possible in the Liouville 
space. The TFD Hilbert space is also specially entangled 
and compressed in order to include asymptotic nonequi- 
librium boundary conditions, in other words, all equilib- 
rium or nonequilibrium statistical mechanics are embed- 
ded directly in the Hilbert space of the system in the 
TFD. 

Let us now consider a system in nonequilibrium steady 
state. A steady state is established by the delicate bal- 
ance between irreversible processes and the driving forces 
produced by the reservoirs. Nonequilibrium steady-state 
systems are ubiquitous in nature and their theoretical de- 
scription has been a challenging fundamental problem for 
many yearsJ^ We define a nonequilibrium steady state as 
an asymptotic state (t — > 0) of a finite quantum system 
which was placed at t = — oo into contact with several dif- 
ferent macroscopic thermal baths or particle reservoirs ^ 
We write the total thermal Hamiltonian in the following 
form 

H(t,r)) = H-H = Hs+H B +eM r lt)'HsB, (8) 

where 77 is a small positive quantity. The system is iso- 
lated from the baths at infinite past, so the system and 
the baths are initially prepared in thermal equilibrium: 

M-oo)} = |Vs> <8 |V>b>, (9) 

where 

\iPs) = -^^expt-^T,)]*) ® \s), (10) 



\i, B ) = [] -4f= E ex P (-£V2T B )|&) ® |6>, (11) 
b V B beB 

and obviously (Hs + "Hb)\iP(— 00)) = 0. What nonequi- 
librium state will develop adiabatically from 00)) 
© if we turn on the system bath coupling? If we follow 
Gell-Mann and Low derivation^ and perform the calcu- 
lations with finite adiabaticity parameter 77, it is possible 
to show that 

(Hs + HB + HsB)\ipr,) = irig-^\ip v ), (12) 

where TLsb is proportional to a coupling strength g and 
\tp v ) = fexp(-i /_ ao dr7<(T,J7))|V'(-oo)). Contrary to 
the ordinary quantum mechanics nonequilibrium wave- 
function does not acquire an infinite phase. Therefore, 
if the system-bath coupling is switched on adiabatically 
(i.e. X] — ► 0), then the r.h.s of eq. (fT2]) becomes zero. It 
is mathematically and physically similar to the obser- 
vation that for the classical nonequilibrium systems the 
slowly switching the interaction on changes the distribu- 
tion but not the eigenvalue of the Liouville operator.^ 
Therefore, if we switch on the coupling between system 
and baths adiabatically slow, the initial equilibrium dis- 
tribution |9]) will evolve to the nonequilibrium steady 
state wavefunction l^) which is an eigenstate of the total 
thermal Hamiltonian with zero eigenvalue. Let us sum- 
marize the formal results: (i) the TFD time-dependent 
Schrodinger equation (j4|) is completely equivalent to the 
Liouville equation for nonequilibrium density matrix, (ii) 
\ip{t)} obtained from (0| contains all information about 
the system necessary to compute any nonequilibrium av- 
erages, (iii) for the system connected to the baths of 
non-interacting particles, the boundary condition that 
the system-environment is decoupled in the initial state 
are equivalent to the requirement that \ip(t)) becomes 
a vacuum for annihilation operators defined via thermal 
Bogoliubov transformation for the bath particles ([7]), (iv) 
the steady state is the state with eigenstate zero of the 
total thermal Hamiltonian of the system. 

Having discussed general theoretical aspects of the ap- 
proach, let us consider an important specific example of 
electron transport through the interacting region. Wc 
begin with the general tunneling Hamiltonian: 

H = H s (al,a s ) + ^ £ b a b a & + $Z {tb s a\a s + h.c), 

beL.R s.beL.R 

(13) 

where a^(ab) creates (destroys) an electron in the state 
b of either the left (L) or the right (R) lead, and a\ and 
a s are creation annihilation operators in the interacting 
region. Here Hs(a\,a s ) represents the interacting re- 
gion (for example, quantum dot, atom or molecule) and 
contains two-particle electron-electron correlations and, 
if necessary, electron-phonon interactions. First, we in- 
troduce the thermal Hamiltonian by formally doubling 
the degrees of freedom H = H — H . Then, we impose 
the reservoir boundary conditions by performing thermal 
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canonical Bogoliubov transformation |J7J| on the particles tonian becomes: 
from the left and the right leads. So the thermal Hamil- 



H = Hs(al,al,a s ,a s ) + 



E 

beL.R 



s.beL.R 



n b {alpl+ alfil) +h.c. 



(14) 



By means of the canonical Bogoliubov transformation 
for reservoir operators, we put all boundary conditions 
inside the thermal Hamiltonian. The thermal Hamilto- 
nian (|14[) depends now on temperatures and chemical 
potentials of the leads. It contains terms as a/3, (3' etc 
which usually appear in the theory of superfluidity (so- 
called "dangerous Bogoliubov diagrams")^ Due to these 
terms the structure of the nonequilibrium steady state 
for the electron transport problem will have some mathe- 
matical similarities with the vacuum in the Hartree-Fock- 
Bogoliubov theory- 
Let us apply Wick theorem to the thermal Hamiltonian 
of the interacting region 

H s (alala s ,a s )=H s +:Hs: (15) 

Here Ti° s is the quadratic part of the thermal Hamiltonian 
and : TLs : is the normal ordered part which contains ir- 
reducible to the quadratic two-particle interactions. The 
normal ordering is performed with respect to nonequilib- 
rium vacuum which is yet to be defined. We partition the 
total thermal Hamiltonian into "simple" part 7io, which 
collects all quadratic terms from the systems, leads and 
system-lead interactions, and the normal ordered part of 
the system thermal Hamiltonian: H = Ha+ ■ H-s '■■ The 
" simple" part of the thermal Hamiltonian H.q can be ex- 
actly diagonalized by some unitary rotation U of annihi- 
lation and creation operators. Then, if we combine this 
unitary rotation U and the thermal Bogoliubov trans- 
formation ||7|), the total transformation from the "bare" 
operators to nonequilibrium quasiparticle operators can 
be written schematically as 

" = W (a!) Wt < ^ 

where (a) are creation (annihilations) operators for 
nonequilibrium quasiparticles. The unitary transforma- 
tion operator W takes into account both, nonequilib- 
rium reservoir boundary conditions (cos 6 b = \f\ — n b ) 
and many particle correlations: 

W = Uexp ( ~ ^beLM °b(aldl - a b a b ) ^ ^ 

The operator exponent in this transformation represents 
the canonical Bogoliubov transformation j7]) performed 
on the operators in the left and the right leads. The total 



thermal Hamiltonian in the nonequilibrium quasiparticle 
basis becomes 

H = ^E k (ala k -a\a k )+: H S :, (18) 

k 

where is the spectrum of nonequilibrium quasiparti- 
cles. The steady state nonequilibrium wavefunction |0) 
for the "simple" part of the thermal Hamiltonian Ho 
(7io|0) = 0) is also a vacuum for nonequilibrium quasi- 
particles a 

a \0) = a\0) = 0. (19) 

Now we know nonequilibrium steady state vacuum and 
corresponding nonequilibrium quasiparticles, we can con- 
struct the nonequilibrium Fock space for our system: 

|0), a+10), 5+10), atfi+10),.... (20) 

We note here that the thermal Hamiltonian, by it- 
self, does not mix original and "tilde" Fock spaces, the 
nonequilibrium mixing is produced by so-called vacuum 
correlation since we introduce the common vacuum Q19p 
for both spaces. The somewhat different idea of defin- 
ing a nonequilibrium Fock space, so called third quanti- 
zation, has been recently proposed by Prosen based on 
quadratic Lindblad master equation. 21 . 

We would like to emphasize here the following impor- 
tant point. Fifty years ago Bogoliubov and Landau in- 
troduced one of the most beautiful and useful concepts in 
theoretical physics, a quasiparticlei 22 i 23 The basic idea of 
quasiparticles is to represent true ground state of inter- 
acting many particle systems as a vacuum with respect to 
some quasiparticle annihilation operators. The quasipar- 
ticle description is, in principle, exact if one defines the 
quasiparticles in terms of the exact eigenstates of many- 
particle systems. Although in practice it can be done 
only approximately by establishing the relation between 
quasiparticles and "bare" particles of the system, for ex- 
ample, by canonical Bogoliubov transformations. Here, 
we demonstrated how one can define quasiparticles for an 
interacting quantum many-particle system in nonequilib- 
rium in such a way that they include both nonequilibrium 
and correlations into their structure. 

With the nonequilibrium Fock space and the ther- 
mal Hamiltonian, which is defined in this Fock space, 
we can readily extend standard methods of electronic 
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structure calculations, such as many-body perturba- 
tion theory, coupled-cluster technique, renormalization 
group, or configuration interaction, to nonequilibrium 
steady state systems. For example, the coupled clus- 
ter equations for nonequilibrium electron transport can 
be easily formally obtained if we make the following 
ansatz for correlated nonequilibrium steady state wave- 
function \ip) — exp(S')lO), where |0) is the vacuum for 
nonequilibrium quasiparticles and S — Sij a l a j + 

J^ijki Sijkia\aja k ai + ... M Here the summation is re- 
stricted to the subspace of the interacting region. Then 
the equations for coefficients S can be obtained by the 
standard methods^ with the use of the thermal Hamil- 
tonian (TTg|) . 

In order to compare our approach with existing the- 
oretical methods, we consider a single impurity con- 
nected to two leads held at different chemical potentials. 
Here we can perform most calculations analytically and 
compare the final expressions with those obtained with 
nonequilibrium Green's functions. We begin with the 
thermal Hamiltonian (|14p where the system comprises 
only one site: 

Hs{at,al,a s ,a s ) = e(a l a — aJa). (21) 

We use Heisenberg equations of motion for all system and 
leads creation/annihilation operators iA = [A,H], where 
TL is the total thermal Hamiltonian. Then we perform 
the Laplace transform f(z) = J Q dt exp(— zt) f (t) for all 
time dependent operators in these equations. Working 
with the Laplace transform it is more convenient to as- 
sume that the system was prepared in equilibrium at 
t = and the nonequilibrium steady state is established 
at t = +oo. In this case we are able to compute di- 
rectly the asymptotic operators by simply multiplying 
the Laplace transformed operator by z and then letting 
z go to zero: lim t ^ +00 f(t) — lim 2 ^o z fi z )- We elim- 
inate all bath operators from the Heisenberg equations 
of motion and the system particle annihilation operator 
becomes 

a(z) — iG(z)a 

+ l G(z)J2^^(VT^7T b p b + 7^,(22) 

b 

where a, /3, and [y in the right side of this equation 
arc time-dependent operators computed at t = 0, a(z) 
is the Laplace transform of the system particle annihi- 
lation operator a(t), and G(z) = [iz — e — E(z)] -1 is 



the system Green's function. S(z), which enters the 
Green's function, is the standard self-energy of the leads 
S(z) = Sl(-z) + T,r(z). Let us now compute, as an ex- 
ample, the steady state electron density on the impurity 

/>oo 

p=\imz dtexp(-zt)(a^(t)a(t)), (23) 
Jo 

where averaging (...) is performed over the vacuum for a, 
(3, and (3 operators. If we substitute the expression for 
a(z) (HU) into (|23"]l. we obtain 

i r + °° 

P=^-J duj\G(iuj)\ 2 X(uj), (24) 

where X(u) = Es=L,_R n ( CJ _ A*b)[( s b(-^) - 
Es(— iw))]. The same result can be obtained by the 
method of nonequilibrium Green's functions through the 
calculation of G-lesser Green's function G< for the sys- 
tem (see, for example^). Likewise, the very similar cal- 
culations can be performed for the current operator (not 
shown here) , and the the standard Landauer formula can 
be recovered. 

In conclusion, based on the formalism of TFD we 
have proposed the concept of nonequilibrium quasipar- 
ticles and nonequilibrium Fock space by representing a 
nonequilibrium steady state of interacting many parti- 
cle systems as a vacuum with respect to some quasi- 
particlc annihilation operators. We demonstrated that 
these nonequilibrium quasiparticles include correlations 
and nonequilibrium into their structure. We developed 
a general theory for the quantum many-particle systems 
connected to several thermal baths or particle reservoirs 
as well as the general theoretical framework to deal with 
electron transport through the interacting regions. We 
have discussed how the standard methods of electronic 
structure calculations such as coupled cluster theory can 
be extended to electron transport problem. On the ex- 
ample of electron transport through the single level impu- 
rity we have demonstrated that the results obtained from 
our approach are consistent with those obtained from the 
nonequilibrium Green's functions. 
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